Program Synthesis of Sparse Algorithms for Wave Function and Energy Prediction in Grid-Based Quantum Simulations

We have recently shown how program synthesis (PS), or the concept of “self-writing code”, can generate novel algorithms that solve the vibrational Schrödinger equation, providing approximations to the allowed wave functions for bound, one-dimensional (1-D) potential energy surfaces (PESs). The resulting algorithms use a grid-based representation of the underlying wave function ψ(x) and PES V(x), providing codes which represent approximations to standard discrete variable representation (DVR) methods. In this Article, we show how this inductive PS strategy can be improved and modified to enable prediction of both vibrational wave functions and energy eigenvalues of representative model PESs (both 1-D and multidimensional). We show that PS can generate algorithms that offer some improvements in energy eigenvalue accuracy over standard DVR schemes; however, we also demonstrate that PS can identify accurate numerical methods that exhibit desirable computational features, such as employing very sparse (tridiagonal) matrices. The resulting PS-generated algorithms are initially developed and tested for 1-D vibrational eigenproblems, before solution of multidimensional problems is demonstrated; we find that our new PS-generated algorithms can reduce calculation times for grid-based eigenvector computation by an order of magnitude or more. More generally, with further development and optimization, we anticipate that PS-generated algorithms based on effective Hamiltonian approximations, such as those proposed here, could be useful in direct simulations of quantum dynamics via wave function propagation and evaluation of molecular electronic structure.


PS function definitions
: Definition of input vectors y and matrices M defined for the n i = 11 input terminals used in PS simulations. Here, V is the PES evaluted on the input coordinate grid. Explicit elements of the matrix are written as M ij , with M ii representing the matrix diagonal elements.

Input index Input vector
Input matrix The elements x i represents the position of the i-th grid-point in the uniform grid. Except where indicated, all operations act in an element-wise manner on all entries in the workspace vector or matrix. We note that I 0 (x) is the Modified Bessel Function (first kind, zero order), as implemented in numpy. Considering all combinations of constants, c = [m, 2, 3, π, 4, L], and functions, we have a total of n r = 134 possible operations at each internal code line. We note that this list of functions is clearly not exhaustive and is somewhat arbitrarily chosen based on trial-and-error investigations; however, the results of the main article demonstrate that this set is sufficient to generate new codes which perform very well in eigenfunction prediction when compared to standard DVR schemes.
It is interesting to ask if the PS-generated algorithms from Table 1 in the main text (which work well across a range of grid-sizes and broadly demonstrate similar convergence to CM-DVR) can be improved upon by seeking separate optimal algorithms for each different gridsize. By removing the constraint that a PS-generated code has to work well for a range of grid-sizes, we might anticipate that generating different algorithms tuned for different grid-sizes might lead to further improvements in accuracy. As such, we performed further PS optimizations but, instead of using a range of grid-sizes in the randomly-generated PES target data, we used single specific grid-sizes. Here, for code size N = 20, we performed 100 PS optimizations for fixed grid-sizes n g = [13, 15, 19, 21, 31] using method E1; beyond fixing the grid-size, the remaining calculation details were the same as noted in sections 2 and 3.
At each grid-size, we then selected the best algorithm from each of the 100 PS optimization runs, and further evaluated the performance at the targeted grid-size by calculating E f for 500 randomly-generated PESs. Figure 1 shows the results of these simulations as a function of grid-size, compared to the CM-DVR method. (working equations given below). It is found that the grid-targeted algorithms generally improve on the codes generated for grid-ranges (Fig. 4), as might be expected; for example, whereas the average RMS fractional errors obtained by the algorithms in Fig. 4 are around 0.5-1.0% smaller than the corresponding CM-DVR results, in the case grid-optimized algorithms, it is found that the errors are decreased further still, reducing the RMS fractional errors relative to converged CM-DVR by up to ∼9%. Of course, the price paid for this small improvement is the requirement of using different algorithms for different grid-sizes, which is not particularly convenient if one is interested in using PS-generated codes in general analysis of quantum molecular vibrational properties. Here, unique algorithms were generated by PS which were specifically optimized to work for a single grid-size; the PS results here (blue circles) show the RMS fractional errors for five different algorithms targeted to work for different grid-sizes. For comparison, the RMS fractional errors for the Colbert-Miller DVR scheme are also shown (red dashed line). Fractional errors were calculated as in main manuscript for the first n e = 3 eigenstates of 500 randomlygenerated PESs; error bars are typically much smaller than the symbol sizes and are not shown for clarity.

Working equations for grid-optimized algorithms
The following working equations were derived from PS simulations optimized for different grid-sizes: • n g = 13 where • n g = 15 where • n g = 19 • n g = 31 where